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Abstract. We investigate the sensitivity of fits of hadron spectra produced in 
heavy ion collisions to the choice of statistical hadronization model. We start 
by giving an overview of statistical model ambiguities, and what they tell us 
about freeze-out dynamics. We then use Montecarlo generated data to determine 
sensitivity to model choice. We fit the statistical hadronization models under 
consideration to RHIC data, and find that a comparison x 2 fits can shed light on 
some presently contentious questions. 
PACS: 12.38.Mh, 25.75.-q, 24.10.Pa 



1. Statistical hadronization models 

There is increasing interest to apply the statistical model [1, 2, 3] in a description 
of both particle yields [4, 5, 6, 7] and spectra [8, 9, 10, 11, 12, 13, 14] produced 
in heavy ion collisions in terms of a statistically hadronizing transversely expanding 
fireball. The fitted parameters, and in particular the temperature, however, have 
varied considerably, ranging from as low as 110 MeV [8, 9, 10, 11, 14] to 140 MeV 
[4, 13] to as high as 160 and 170 MeV [5, 6, 7, 12]. This report examines these 
differences in terms of the different possible approaches to statistical hadronization, 
and studies data sensitivity to model choice. 

The most commonly used prescription for statistical emission of hadrons is the 
truncated Cooper- Frye formula [15]. 

( E ^) cf = / P"d%J( P W,T, A)0(p"d%) + ( E ^) Res (!) 

Where p^ is the particle's four-momentum, is the system's transverse flow, T 
is the temperature, A is the normalization and/or chemical potential, f(E,T,\) is 
the statistical distribution of the emitted particles in terms of energy and conserved 
quantum numbers, S M describes the hadronization geometry and 8(p f ^d 3 'E f j,) is a step 
function, which corrects for inward emission [16, 17]. The Cooper- Frye formula is 
appealing because it's the most general covariant way of describing statistical particle 
emission. However, its form is compatible with a variety of freeze-out scenarios. 

1.1. Resonances and normalization 

As Fig. 1 shows, the resonance contribution (E^^)^ es is a non-negligible and strongly 
temperature dependent admixture to the total multiplicity. Its role within particle 
spectra, however, is a matter of controversy. If an interacting hadron gas phase re- 
equilibrates resonance decay products [18], it might be more realistic to not include 
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Figure 1. Left: Fraction of particles produced by resonances as a function of 
Temperature, Right: dNtot / d(y , rf) arising from freeze-out of a Bjorken fluid, 
before (open circles) and after (solid circles) resonance decays. The normalization 
is arbitrary, chosen to coincide with PHOBOS [28] data (triangles) 



resonances when analysing spectra [8, 9, 10, 11]. Direct detection of resonances 
[19, 20, 21, 22] through invariant mass reconstruction, however, suggests that at least 
in some cases their feed-down has to be taken into account. These resonances will have 
a different distribution than those produced directly, and their inclusion can shift the 
temperature by as much as 20 MeV [23]. 

If an incoherent system is considered, one can factor out the decay matrix elements 
(the angular distribution of decay products and lifetime) , and arrive at the distribution 
of "daughter" particles (N with momentum p) from "parent particles" (Nr with 
momentum pn) by simply integrating over the available phase space 

In general, this expression can get very complicated, and Montecarlo integration [24] 
becomes necessary. For most cases considered here, where there is one feed-down and 
two or three body decays, eq. 2 can be integrated semi-analitically [12, 13, 25]. 

While including resonances in this way is computationally demanding [25] , it can 
be implemented with no extra degrees of freedom. It is therefore straightforward to 
check the effect of resonance inclusion on \ 2 ■ A large increase in statistical significance 
would be a strong indication that resonance decay products do in fact emerge from the 
system with little rescattering. Otherwise, rescattering has to be taken into account. 

The handling of resonance decays is connected with the way particle spectra are 
normalized. Several ways to normalize spectra are possible, and the fitted temperature 
could change significantly depending on the chosen method, since temperature 
becomes correlated to a different set of parameters. 

The simplest approach is to treat normalization of each particle as a free 
parameter [8, 9, 10, 11]. If the post-hadronization interacting hadron gas phase is 
long enough to allow inelastic interactions to alter the chemical composition of the 
fireball, this may be the correct approach. In this case, temperature will be correlated 
to flow but not to normalization. 
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The arbitrary normalization means there is no way to include feed-down from 
resonances without introducing many more parameters (the resonance normalizations) 
into the fit. (If rescattering is significant, however, these parameters would be 
physically justified). 

Alternatively, a common normalization volume, together with the introduction of 
flavor chemical potentials, can be used [12, 13]. This approach, justified by the success 
in handling short-lived resonances through thermal models [7], is consistent with a 
scenario in which post-hadronization dynamics does not change particle abundances or 
distributions significantly. It does not require extra degrees of freedom for resonances 
to be included. 

However, feed-down from resonances means the fitted temperature will be 
strongly correlated with chemical potentials and reaction volume, in addition to 
transverse flow: Most short-lived resonances have the same quark content as the 
lighter daughter particles, making their yield relative to the lighter particle dependent 
on temperature only. Hence, temperature will influence both absolute normalization 
(correlating it with volume, chemical potentials) and the slope (correlating it with 
flow). 

If chemical potentials are used, the eventuality of chemical non-equilibrium 
raises the possibility that, additionally to chemical potentials, the flavor phase space 
occupancy parameter 7 needs to be used for cither strange [6, 12] or light [4, 26] flavors. 
If these parameters are used, they will correlate temperature with chemical potential 
even more strongly, since, just like temperature, 7 affects particles and antiparticles 
abundance in the same way. 

1.2. Freeze-out geometry and flow profile 

Several choices of the freeze-out geometry are possible. This choice is important 
both for getting a full understanding of how freczc-out happens [27] and because 
it correlates with the shape of particle spectra, and hence on parameters such as 
transverse flow and temperature. 

Measured RHIC rapidity distributions [28, 29] indicate that around midrapidity 
Bjorkcn boost-invariance [30] holds. This, together with cylindrical symmetry 
(appropriate for the most central collisions) constrains rf 3 S M to be of the form [25, 27] 

^ = [f / cosh(y L ),rcos(6»),rsin(6»),t / sinh(y L )] (3) 

dt f 

p^d'En = (m T cosh(y - y L ) - -^p T cos(# - 4>))rdrdy (4) 

or 

where y, <f> are the particle's rapidity and momentum directions, while yL,Q are the 
longitudinal rapidity and emission angle of the freeze-out hypersurface. t/ is the 
freeze-out time in a longitudinally co-moving frame. 

It is clear that several freeze-out models can be produced with different choices 
of dtf/dr, giving quantitatively different results. It is also clear that the choice of 
dtf /dr might influence fitted parameters such as temperature, transverse velocity and 
normalization. Recently published fits use a variety of freeze-out geometries: 

• [8, 9, 10, 11] assume dtf /dr — 0, i.e. emission is independent of radius. Emission 
happens in the same lab momentum time. 

• [13] fits dtf/dr as a free parameter, with the result being ~ 1. 

This indicates explosive freeze-out, driven by negative vacuum pressure [31] 
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Figure 2. Results of fits to Montecarlo generated data samples 1 and II (-p- = 0): 
X 2 profiles for temperature (left) and fitted dtf/dr (right). The models used in 
the fits are described in sections 1.1 and 1.2 

Full fit includes fitted chemical potentials and 7 9 , s , resonances and fitted dtf/dr. 
Last profile shows effect of fitting sample II using Eq. 6 (Sample I Full fit) 



• In [12] freeze-out happens at the same proper time r throughout the fireball, 
(g? 3 £ m oc Freeze-out might have this form ("Hubble" freeze-out) if dynamical 
effects are negligible during hadronization [27]. 

An additional physical uncertainty, which is found to be correlated to flow and 
dtf/dr in fits, is the flow profile. Hydrodynamics requires that each part of the fireball 
volume will in general have a different density and transverse expansion rate. For this 
reason, the integral over <i 3 £ will in general span a range of flows, weighted by density 
function p(r). 

E 1T = I rdr ( E ~ PT^r)f{T, y T (r), A)p(r) (5) 
d°p J dr 

Just like the form of <i 3 £, the flow profile depends on the hadronization conditions. For 
instance, assuming that freeze-out happens at a critical energy density yields yr oc r 
[32] , while assuming freeze-out happening at approximately the same time in the lab 
frame leads to a power distribution for the transverse expansion velocity vt oc r a [33] 
One can approximate the flow profile with one "average" flow and assume [10, 13] 

d,N dtr 
E^-^(E-p T ^l)f(T,<y T >,X) (6) 
d°p dr 

However, this approach might result in a systematic shift of fitted parameters such as 
flow and dtf/dr 

2. Sensitivity to model choice 

The ambiguities presented in sections 1.1 and 1.2 mean that it is important to study 
their effect on the statistical model's fitted parameters. One way to do this is to use 
a Montecarlo to generate data according to a particular freeze-out model, and to see 
what happens if the "wrong" model is used to perform the fit. We have written a 
MonteCarlo program which can be used for this purpose. While a detailed description 




Figure 3. Results of fits to Montecarlo generated data samples III (-p- = 1): 
X 2 profiles for temperature (left) and fitted dtf/dr (right). The models used in 
the fits are described in sections 1.1 and 1.2 

Full fit includes fitted chemical potentials and 7 9lS , resonances and fitted dtf/dr. 



of the algorithm is outside the scope of this write-up [34], we shall provide a short 
summary of its features here. An acceptance/rejection algorithm is used to generate 
particles in a statistical distribution in the volume element's rest frame. The accepted 
particles are then Lorentz transformed to the lab frame. (Any flow and density profile, 
as well as any freeze-out surface can be accomodated) . Resonance decays are handled 
through Eq. 2, using the MAMBO algorithm [24] to generate points in phase space. 
Output can be used to generate spectra or fed into a microscopic model such as 
uRQMD [35]. 

The Montecarlo output was used to produce the data points in Fig. 1. It can be 
seen that a boost-invariant statistical hadronization can explain the global properties 
of the system such as dNtot/drj. It can also be verified that the role of resonances is 
absolutely crucial. 

We proceeded to generate three datasets of particles. Each data set had a 
temperature of 140 MeV, a maximum transverse flow of 0.55, and out of equilibrium 
chemistry (7^ = 1.4, "f s Hq — 0.8). Generated particles include 7r, K + , K~ 7 p,p, A, E, Q 
and their resonances. The three samples differ in their choice of freeze-out geometry 
(specifically dtf/dr) and flow profile: 

Sample I dtf/dr = and no flow profile, as fitted in [10] 

Sample II dtf/dr = and a quadratic flow profile, as fitted in [8, 9, 11] 

Sample III dtf/dr = 1, the boost-invariant analogue of [13]. 

We have fitted the three samples to a variety of models, producing x 2 /DoF profiles for 
freeze-out temperature and the dtf/dr parameter. Fig. 3 shows the profiles resulting 
in the fit to sample III while samples I and II are shown in Fig. 2. A full chemistry 
model with resonances (solid black line) seems to be equivalent (as far as the position 
of the x 2 minimum and the value of x 2 / 'DoF) to a fit in which normalization is 
particle-specific (dashed red line). However, the chemical potential fit has greater 
statistical significance since considerably more degrees of freedom are required for 
arbitrary normalization. 
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Figure 4. Fit to RHIC data 

profiles for temperature (left) and fitted dtf/dr (right). The models used in the 
fits are described in sections 1.1 and 1.2 

Full fit includes fitted chemical potentials and 7<j,s, resonances and fitted dtf/dr. 



If chemical potentials are included resonances become essential since a fit 
with chemical potentials but no resonances (dotted blue line) loses all statistical 
significance. Similarly, the physical presence of non-equilibrium (7 9 , s ^ 1) means 
chemical potentials have to include the non-equilibrium parameter for the fit to be 
meaningful (black dot-dashed line) . The freeze-out geometry does not seem to impact 
the temperature minimum that much. However, the correct freeze-out geometry can 
be picked out by a comparison of fits to different models by choosing the model with 
the lower y^/DoF. Moreover, data sensitivity to temperature is strongly affected 
by freeze-out geometry: Comparing the temperature profiles for different choices of 
dtf/dr (Figs. 2 and 3) it is apparent that the temperature x 2 minimum is more 
definite in the dtf/dr = case. In the case of explosive freeze-out, the correlation 
between temperature and other parameters in the fit (notably flow) increases, resulting 
in a shallow \ 2 increase at larger than minimum temperatures. 

Finally it should be noted that flow profile, freeze-out geometry and temperature 
appear to be strongly correlated. If data sample II is fitted with a distribution with 
no flow profile (such as Sample I) there is a non-negligible shift in both the fitted 
temperature and dtf/dr, and a small rise in \ 2 /DoF (Fig. 2, brown dot-dashed line). 

3. Fit to RHIC data (y/s^N = 130GeV) 

Finally, we have perfomed a fit to the available RHIC data. The data sample we used 
is the same as the one used for the Montecarlo, but, since the STAR and PHENIX 
spectra had different trigger requirements (notably centrality) and acceptance regions, 
we have used two different system volumes, one for STAR particles and another one 
for PHENIX. 

The results are similar to the Monte Carlo data in many ways. The \ 2 /DoF was 
only slightly larger. 

A fit with particle-specific normalization gives a very similar \ 2 and fitted temperature 
to a fit including chemistry and resonances. If chemistry is included in particle spectra 
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analysis, resonances and non-equilibrium are essential. The fit to freeze-out geometry 
weakly points to dtf/dr = 1, a picture that is supported by the temperature x 2 
profile, virtually identical to data sample III. However, we can not claim our study to 
be complete in this respect, since we have not yet investigated the effect of including 
flow profile in the models. As Montecarlo simulations have shown, the conclusion can 
differ once these are taken into account. 

Fig. 5 shows the hyperon and pion spectra from the global fit of Fig. 4. A 
comparison of the fits on the left panel confirms that a model with no chemistry is 
about as good at fitting particle spectra as a model with resonances and chemical 
potentials. However, the second fit also has predictive power: We have used the fitted 
parameters to predict the tot spectrum for the £*. Unsurprisingly, we found that 
the S* should have roughly the same slope as the EEs, but its total multiplicity should 
be about three times as big. We therefore suggest that a greater sample of spectra, 
in particular more spectra of heavy resonances taken within the same centrality bin 
as light particles (to make sure both flow and emission volume are the same for each 
particle) would help in establishing whether chemical potentials are a good way to 
normalize hadron spectra or not. 

The only spectrum which presents a significant systematic deviation for most 
models is the tt~ . As Fig. 5 (right panel) shows, most models fail to catch the upward 
dip of the low momentum pions, and indeed simple transverse expansion predicts the 
reverse trend [8, 9, 10, 11]. Including resonances, and allowing for -f q > 1 helps (the 
latter is equivalent to postulating a pion "chemical potential" [36]) However, to fully 
account for the lowest momentum pions, even addition of resonances and 7 9 > 1 are 
not enough. One has to add a power-law component to the pion spectrum 

p dN = ( dN\ A 

^d^p \ <Pp) c f (PT+PT0) a [> 

This contribution (roughly 6% of the total pion yield in the best fit) also accounts 
for the highest pt pions. Such a parametrization has been justified [37, 38] and 
successfully used [39] before. 

In conclusion, we have presented a comparison of different statistical models 
currently used to fit spectra in heavy ion collisions. We have described how these 
different models arise from different freeze-out scenarios of a system hadronizing from 
a thermalized quark gluon plasma. We have used MonteCarlo simulated data to study 
the sensitivity to model choice of presently available experimental data, and have 
evaluated different models ability to fit presently existing RHIC data. While data 
slightly favors a chemical non-equilibrium explosive freeze-out, there is not enough 
evidence to make a definitive conclusion about this issue. We hope the forthcoming 
200 GeV results will clarify this further. 
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Figure 5. Left: S and Q tut distributions, together with the global fit (from 
Fig. 4, solid black and dashed red). The model with chemical potentials (solid 
black) was used to predict the E* mj- distribution using the Montecarlo. 
Right: PHENIX ir~ pt distribution, with the global fit from Fig 4. As can 
be seen, both resonances and ■yq help, but are not sufficient to explain the pion 
distribution fully. 



[2] I. Pomeranchuk, Proc. USSR Academy of Sciences (in Russian) 43 (1951) 889. 

[3] R. Hagedorn, Suppl. Nuovo Cimento 2 (1965) 147. 

[4] J. Rafelski and J. Letessier, arXiv:nucl-th/0209084. 

[5] P. Braun-Munzinger, I. Heppe and J. Stachel, Phys. Lett. B 465, 15 (1999). 

[6] S. V. Akkelin, et al, Nucl. Phys. A 710, 439 (2002) [arXiv:nucl-th/0111050]. 

[7] D. Magestro, J. Phys. G 28, 1745 (2002) [arXiv:hep-ph/0112178] . 

[8] J. M. Burward-Hoy for the PHENIX collaboration, QM 2002 and proceedings 

[9] M. Van Leeuwen for the NA49 collaboration, QM 2002 and proceedings 

[10] Ladislav Sandor for the NA57 collaboration, these proceedings 

[11] J Castillo for the STAR collaboration, these proceedings 

[12] W. Broniowski and W. Florkowski, Phys. Rev. Lett. 87, 272302 (2001) [arXiv:nucl-th/0106050]. 

[13] G. Torrieri and J. Rafelski, New J. Phys. 3, 12 (2001) [arXiv:hep-ph/0012102]. 

[14] K. A. Bugaev, M. Gazdzicki and M. I. Gorenstein, arXiv:hep-ph/0211337. 

[15] F. Cooper and G. Frye, Phys. Rev. D10 (1974) 186 

[16] K. A. Bugaev, Nucl. Phys. A 606, 559 (1996) [arXiv:nucl-th/9906047]. 

[17] C. Anderlik et al, Phys. Rev. C 59, 3309 (1999) [arXiv:nucl-th/9806004]. 

[18] M. Bleicher and J. Aichelin, Phys. Lett. B 530 (2002) 81 [arXiv:hep-ph/0201123]. 

[19] P. Fachini [STAR Collaboration], arXiv:nucl-cx/0211001. 

[20] G. Van Buren for the STAR collaboration, QM 2002 and proceedings 

[21] C. Markert [STAR Collaboration], J. Phys. G 28, 1753 (2002). 

[22] V. Friese [NA49 Collaboration], Nucl. Phys. A 698 (2002) 487. 

[23] G. Torrieri and J. Rafelski, J. Phys. G 28, 1911 (2002) [arXiv:hep-ph/0112195]. 

[24] R. Kleiss and W. J. Stirling, Nucl. Phys. B 385, 413 (1992). 

[25] E. Schnedermann et. al Phys. Rev. C 48, 2462 (1993) [arXiv:nucl-th/9307020]. 

[26] J. Letessier and J. Rafelski, Int. J. Mod. Phys. E 9, 107 (2000) [arXiv:nucl-th/0003014], 

[27] G. Torrieri and J. Rafelski, arXiv:nucl-th/0212091. 

[28] A. Olszewski et al. [PHOBOS Collaboration], J. Phys. G 28, 1801 (2002). 

[29] P. Staszel et al, Acta Phys. Polon. B 33, 1387 (2002). 

[30] J. D. Bjorken, Phys. Rev. D 27, 140 (1983). 

[31] J. Rafelski and J. Letessier, Phys. Rev. Lett. 85, 4695 (2000) [arXiv:hep-ph/0006200]. 

[32] D. Teaney, J. Lauret and E. V. Shuryak, arXiv:nucl-th/0110037. 

[33] E. Schnedermann and U. W. Heinz, Phys. Rev. C 47, 1738 (1993). 

[34] QM2002 poster presentation, Abstract at http://alice-france.in2p3.fr/qm2002/ 

[35] S. A. Bass et al, Prog. Part. Nucl. Phys. 41, 225 (1998) [arXiv:nucl-th/9803035]. 

[36] B. Tomasik, et. al, Heavy Ion Phys. 17, 105 (2003) [arXiv:nucl-th/9907096]. 

[37] R.D. Fields and R.P. Feynman, Phys. Rev. D 15, 2590 (1977); 

[38] R. Hagedorn, Riv. Nuovo Cim. 6N10, 1 (1984). 

[39] T. Peitzmann, arXiv:nucl-th/0303046. 



